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t.O  INTRODUCTION 


Zienkiewicz  (Ref.  1)  and  a  number  of  other  investigators  (Refs.  2 
through  5)  have  noted  that  the  finite  element. method  (FEM)  is  suffi¬ 
ciently  general  to  be  applied  to  a  wide  class  of  problems,  including 
fluid  mechanics.  Within  recent  years,  the  FEM  has  been  applied  to  the 
study  of  certain  fluid  mechanics  problems.  For  example.  Baker  and 
Manhardt  (Ref,  2)  have  applied  the  FEM  to  prediction  of  low-speed 
flows,  and  Chan,  et  al.  (Ref.  5)  have  used  the  FEM  in  the  solution  of 
the  small  disturbance  equation  for  transonic  flow. 

It  is  the  purpose  of  this  report  to  present  results  of  a  study  of  the 
application  of  the  FEM  to  the  compressible  flow  equations  in  the  non¬ 
conservative  Ekilerian  form.  The  flow  over  an  airfoil  in  the  transonic 
regime  is  a  problem  of  current  interest  and  was  chosen  to  demonstrate 
the  power  of  the  method.  The  choice  of  the  Eulerian  equations  as  the 
governing  equations,  rather  than  thin  airfoil  theory,  was  made  to 
alleviate  the  restriction  to  thin  airfoils  and  to  allow  for  extension  to 
other  fluid  mechanics  problems.  Although  solution  to  these  governing 
equations  requires  a  longer  convergence  time  and  additional  machine 
storage  than  does  solution  of  the  one  equation  of  small  disturbance 
theory,  this  generality  will  be  of  future  advantage.  Choice  of  the  non- 
conservative  form  of  the  governing  equations  was  made  for  computational 
convenience  as  the  FEM  algorithm  for  the  conservative  form  is  more 
complex  and  requires  additional  calculations. 

Numerical  results  are  presented  in  this  study  for  subsonic  and 
transonic  inviscid  flow  over  two-dimensional  nonlifting  airfoils.  The 
results  of  the  predictions  are  compared  to  experimental  results  and  to 
other  numerical  solutions. 


2.0  ANALYSIS 


2.1  GENERAL  DISCUSSION 

The  Galerkin  method  of  weighted  residuals  was  used  with  the  un¬ 
steady  Euler  equations  to  obtain  a  system  of  ordinary  differential 
equations  in  time.  The  elements  used  in  the  study  were  parametric, 
bilinear  quadrilateral  elements  (Ref.  1).  During  this  study  various 
diffictilties  were  encountered  with  the  straighforward  application  of 
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finite  elements  to  these  equations.  One  of  the  first  problems  encountered 
was  the  irregular  behavior  of  the  solution  near  the  leading  and  trailing 
edges.  This  was  caused  by  the  discontinuity  in  the  derivatives  of  the 
boundary.  The  following  methods  of  improvement  were  evaluated: 
altering  the  geometry  to  smooth  the  discontinuity,  enforcing  stagnation 
conditions  at  these  points,  'splitting  the  discontinuous  boundary  condi¬ 
tions  and  applying  them  separately  to  corresponding  elements,  and 
altering  the  functional  form  of  the  interpolation  function. 

Another  problem  encountered  was  instability  in  the  time  integration. 
Initially,  an  artificial  viscosity  term  was  used  in  an  attempt  to  remove 
the  instability.  Steady-state  solutions  were  obtained,  but  a  small  error 
was  present  because  of  the  viscosity  term.  The  next  modification  was 
the  periodic  altering  of  the  solution  by  use  of  Bernoulli's  equation. 
Although  the  true  transient  behavior  was  lost,  the  steady-state  solution 
was  obtained  without  using  artificial  viscosity  for  subsonic  cases.  In 
the  transonic  cases,  while  viscosity  was  still  required,  a  steady-state 
solution  could  be  obtained  using  a  very  small  value  for  the  viscosity. 

In  addition  to  the  need  for  the  artificial  viscosity  term  in  the  transonic 
calculations,  a  fine  grid  network  in  the  region  of  the  shock  wave  was 
required. 

2.2  BASIC  EQUATIONS 


The  unsteady  Euler  equations  are 


* 
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The  variables  are  nondimensionalized  by 
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Using  the  isentropic  relation  in  nondimens ional  form  P  ~  p'^  ^  Eq,  (2) 
becomes 

Pi  +  P“x  +  ''f'y  =  0 


y  p, 


UU^  +  VU„  + 
*  / 


P  U 

r  oo  ( 


0 


y  p, 


VVy  + 


jpy-^Py  =  0 


<3) 


The  density  p  is  replaced  as  a  dependent  variable  with  T  = 
yielding 


+  (y  -  DTu,  +  uT^  +  (y  -  DTVy  +  vT^ 

y  Poc 


=  0 


Uj  +  UU^  +  VUy  + 


T  =  0 


(y  -  ‘ip™” 

V  P 

^  oo 


2  X 


Vj  +  uVj^  -I-  vVy  t 


T.  =  0 


(4a)- 

(4b) 

(4c) 


These  are  the  equations  which  are  used  in  this  report.  In  cases 
where  a  shock  occurred  an  artificial  viscous  term  was  needed  to  obtain 
stability.  The  basic  equations  become 


Tj  +  (y  -  l)Tuj^  +  uT^  +  (y  _  1)T  +  vT  =  0 


J 


+  UUj^  +  VUy  + 


(y  -  DM 


2  * 


aV  u 


Vj  +  uv^  +  vv^, 


1 

<y  - 


(5) 


2.3  FORMULATION  OF  THE  ELEMENTS 


The  domain  is  divided  into  small  quadrilateral  elements,  e.  g.  ,  see 
Figs,  1  and  2.  Two  numbering  systems  are  used.  One  is  a  global  sys¬ 
tem  which  assigns  a  sequential  node  number  to  each  node  in  the  domain. 
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The  method  used  to  assign  the  global  system  will  not  affect  the  results 
but  is  a  significant  factor  in  computational  efficiency.  The  other  is  an 
elemental  system  where  the  four  nodes  of  each  element  are  sequen¬ 
tially  assigned  within  each  element.  The  order  used  by  the  authors 
was  counterclockwise  starting  with  the  lower  left  node.  Throughout 
this  report  the  numbering  system  used  will  be  evident  by  the  presence 
of  an  upper  limit  of  4  or  N.  Within  each  element,  the  variables  are 
approximated  as  a  function  of  the  values  at  the  four  nodes  (corners) 
in  the  following  parametric  form  (Ref.  1): 

4 

4 

1=  I 

T  - 

4 

X  =  ,2  Xjfji  (Irj) 

where  (-1  <  ?  <  1  and  -1  <  p  S.  1) 

=  (1  +  m  +  tj)/4 

=  (1  -  ad  +  r,)/4 

0®  =  (1  -  ad  ~  7)/4 

=  d  +  ad  -  a/4  <7) 

The  dependent  variables  can  therefore  be  approximated  over  the  entire 
domain  as  a  linear  combination  of  values  given  at  the  nodes. 

Partial  derivatives  within  an  element  may  be  obtained  from  Eqs.  (6) 
and  (7)  as  follows.  Consider  for  example 

“x  ^  +  “r;'7x 
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where 

vf'J 

J  -  -f*,  - 

giving 


*  ,?i  J.  •' 


4  k' 


k  ^rj  ^7?  I 


k=l 


-  J, k!, 


(8) 


Thus  the  partial  derivatives  are  a  linear  combination  of  the  nodal 
values.  It  is  noted  that  while  the  dependent  variables  are  continuous 
between  elements,  the  partial  derivatives  are  discontinuous  at  the 
boundaries. 

The  time  derivatives  are  treated  as  follows: 


4 

I 

i=l 


=  S  VjQ^ 
i=l 

4  . 

=  s,  T; 

i=i 


(9) 


2.4  GALERKIN  METHOD  OF  WEIGHTED  RESIDUALS 

Substituting  the  approximations  given  by  Eq,  (6)  into  Eq.  (4)  for 
each  element  in  the  domain  gives  equations  of  the  following  form: 

!«,T,  .  F,(T  .j.vj) 

1=1 

S  (ZjUi  =  u.,  v) 


S  ajVi  =  Fgdj,  u..  V  )  ,  j  =  l,  2,  3,  ...N 


(10) 
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where  ori  and  Fj  are  implicit  functions  of  the  independent  variables  x 
and  y  over  the  domain.  The  above  consists  of  three  continuum  sets  of 
equations  with  a  finite  (three  times  the  number  of  node  points)  number 
of  unknowns.  Since  all  of  the  implied  equations  cannot  be  satisfied,  a 
selected  set  of  weighted  integrals  of  these  equations  was  chosen: 


N 


/ 

R 

s 

i=l 

ajTidA  = 

=  l/Ai- 

Uj.  V.)0^dA 

(11a) 

N 

2 

i=l 

QjUj  dA  = 

Uj.  Vj)9i|^dA 

(11b) 

N 

1 

1=1 

ajV.  dA  = 

Uj,  dA 

(11c) 

where  j,k  =  1,2,  . . . ,  N 

The  are  the  basis  functions  formed  from  the  elemental  f^'s. 
There  is  a  basis  function  associated  with  each  global  node,  such  that 
the  function  is  unity  at  this  node  and  zero  at  all  other  nodes.  The  basis 
function  is  a  composite  of  the  elemental  as  follows.  For  a  point 
in  the  domain  that  lies  in  an  element  containing  the  k*^  node,  the  value 
of  is  defined  as  the  value  of  the  n  in  that  element  which  has  a  unit 
value  at  the  node.  The  term  is  defined  as  zero  at  all  other  points. 

In  actual  calculations  the  integration  represented  in  Eq,  (11)  is  done 
separately  for  each  element  and  then  assembled  to  form  the  equations. 
When  the  integration  is  performed  over  a  single  element  only  four  basis 
functions  will  give  nonzero  results.  Within  this  element  these  four 
basis  functions  are  identical  to  the  four  interpolation  functions 
and  The  elemental  system  is  therefore  defined  as 


/p  uj,  v^ni  dA 

m 

/jj  “j’ 

tn 

4  ^3^4’ “j’ 


i  =  1.  2,  3,  4 


(12) 


where  bj  and  Gj  are  implicit  functions  of  the  independent  variables  x 
and  y  over  the  elemental  domain.  The  above  is  a  system  of  twelve 
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equations  in  twelve  unknown  derivatives.  A  more  detailed  analysis  of 
Eq.  (12)  is  shown  in  Appendix  A,  During  assembly,  the  elemental 
systems  are  combined  to  form  Eq.  (11).  This  is  done  by  summing 
together  corresponding  elemental  equations  such  that  the  fl's  in 
Eq.  (12)  form  the  0^'®  same  time,  terms  that 

correspond  to  the  same  unknown  are  combined,  enforcing  continuity 
between  elements.  This  assembly  process  converts  the  individual 
elemental  system  of  Eq.  (12)  into  the  global  system  given  by  Eq,  (11). 

2,5  BOUNDARY  CONDITIONS 

The  examples  used  in  this  report  are  flow  around  nonlifting  two- 
dimensional  airfoils.  The  boundary  conditions  for  these  problems  are 
as  follows.  Consider  Fig.  1  with  the  top  boundary  a  far-removed  wall 
and  the  upstream  and  downstream  boundaries  being  sufficiently  re¬ 
moved  to  assume  free -stream  conditions  at  these  boundaries.  The 
boundary  conditions  are 

u  =  1,  T  =  1,  V  =  0  at  the  upstream  and  dovwistream  boundaries 

V  =  0  at  the  upper  boundary  and  along  the  centerline 

V  =  u  tan  d  along  the  airfoil  (13) 

The  boundary  conditions  are  applied  in  time  derivative  form;  that 
is,  Eq,  (13)  is  differentiated  with  respect  to  time,  obtaining 


'  =  T=  0 

at  the  upstream  and  downstream  boundaries 

(14a) 

v  =  0 

at  the  upper  boundary  along  the  centerline 

(14b) 

=  ii  tan  B 

along  the  airfoil 

(14c) 

If  Eq.  (13)  is  satisfied  for  the  initial  conditions,  then  Eq.  (14)  will 
enforce  Eq.  (13)  for  all  time. 

When  Eq.  (14a)  is  applied  to  a  node,  the  corresponding  equation 
from  Eq.  (lib)  is  removed  from  the  system.  Likewise  when  Eq.  (14b) 
is  applied  at  a  node,  a  corresponding  equation  from  Eq.  (11c)  is  re¬ 
moved,  In  the  case  of  Eq.  (14c),  the  equations  in  Eqs,  (11b)  and  (11c) 
are  replaced  by  Eq.  (14c)  and  a  corresponding  linear  combination  from 
Eqs.  (11b)  and  (11c).  The  linear  combination  is  based  on  the  geometry 
at  the  point  and  is  of  the  form 

(Eq.  lib)  +  Ian  $  (Eq.  11c) 
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The  points  at  the  leading  and  trailing  edges  require  special  hand¬ 
ling  in  order  to  treat  the  discontinuity  of  the  slopes  at  these  points.  At 
these  points  the  boundary  conditions  are  applied  to  the  elemental  equa¬ 
tions  that  form  Eq,  (11)  before  assembling.  In  this  manner,  a  different 
value  for  the  boundary  condition  may  be  applied  to  each  side  of  the  same 
point.  Because  of  the  discontinuity  in  slope,  a  discontinuity  in  u  and  v 
also  occurs  at  these  points.  "■  The  assembly  process  is  altered  at  these 
points  to  enforce  continuity  of  T  and  W.  The  above  process  is  easily 
accomplished  as  follows.  The  columns  corresponding  to  u  and  v  at  the 
singular  point  are  combined  into  a  single  W  column  by 

W  column  =  cos  5  (u  column)  +  sin  0  (v  column) 

The  rows  for  u  and  v  at  the  singular  point  are  combined  into  a  single 
W  row  by 


W  row  =  (ii  row)  +  tan  6  (v  row) 

The  previous  operations  result  in  a  new  elemental  system  of  equa¬ 
tions  with  one  less  equation  and  unknown.  In  order  to  maintain  a 
simple  assembly  process  a  dummy  equation  and  unknown  are  added  to 
the  element  system.  This  transformed  elemental  system  is  then  assem¬ 
bled  with  the  rest  of  the  elemental  systems.  The  regular  boundary  con¬ 
ditions  are  then  applied  at  the  other  boundary  nodes.  When  the  global 
system  is  solved,  the  unknowns  at  the  singular  points  correspond  to  T, 

W ,  and  a  dummy  unknown  that  is  ignored.  A  value  for  W  is  given  at  the 
starting  time  and  W  is  integrated  in  time  for  the  singular  points.  The 
velocity  components  at  the  point  are  obtained  by  using  the  value  of  tan  0 
corresponding  to  the  element  under  consideration.  The  values  of  u  and 
V  will  be  discontinuous  at  the  singular  point  and  along  the  line  segment 
connecting  the  two  elements.  They  are  continuous  at  the  other  node 
common  to  both  elements . 


‘The  actual  solution  to  the  equations  gives  u=v~Q  at  these  points 
with  singular  derivatives.  Imposing  u=v=0  as  boundary  conditions  re¬ 
quired  either  a  large  number  of  extra  node  points  in  the  vicinity  of  the 
discontinuity  or  the  inclusion  of  a  special  singular  function  for  the  fit 
at  these  points.  The  method  described  in  the  text  gave  comparable 
results  with  less  computational  effort. 
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2.6  CALCULATION  OF  TERMS  IN  THE  ELEMENTAL  SYSTEM 

OF  EQUATIONS 

The  elemental  system  of  Eq.  (11)  can  be  obtained  in  a  straight¬ 
forward  but  tedious  manner.  For  those  readers  not  familiar  with  the 
mechanics  of  the  finite  element  method,  the  steps  for  obtaining  the 
elemental  equations  for  axial  momentum  [Eq.  (4b)]  will  be  outlined. 
The  other  two  equation  forms  can  be  obtained  in  a  similar  manner. 

The  axial  momentum  Eq.  {4b)  when  written  in  detail  for  an  element 
is 

i  i.n'  .  i  u,  a'  x  u  ni  -  i  v,  n‘  i  .  ni  -  k  x  t,  a', 

1=  1  i=  I  J=  1  i=:  1  j=  1  '  1=  1  V  i  3 1 


Multiplying  by  and  integrating  over  the  element  results  in 


4  4 


i=i  j-i  '  '  'r.  " 

-  I  i  vju  /  -  k  i  tj  n'o^^dA 

^  Rm  ^  i=l  ‘  R..  * 


(16) 


Note  that  while  the  right-hand  side  is  a  function  of  the  dependent  vari¬ 
ables,  the  integrals  are  not;  they  can  be  evaluated  on  a  one-time  basis. 
This  is  true  for  the  last  term  because  T  is  used  as  a  dependent  variable 
instead  of  p.  The  above  integrals  are  integrated  in  the  parametric 
coordinate  system  (?  ,  r) ).  The  integrands  are  polynomials  of  degree 
three  or  less  in  either  f  or  rj  .  Consider  for  example  the  first  term 
on  the  right-hand  side  for  i=j=k=l; 


//  «‘ti^nMxdy  = /;  ;,^)n’dxdy 

R|r  Rm 


=  /■"  ^  n f  y,,  -  y^) ' J1  n  ^ J  dfdr? 

=  +  ad  +  t,)/4][(1  +  r,)y^  -  (1  + 


(17) 


[(1  +  ad  +  r,)/4U^d7, 


(1  +  -  Tf)- 

16 


-1  -1 


td  +  Tj)y^  -  (1  +  ay^  d^drj 


13 


AEDC-TR-76-86 


where 


Yri  =  -  O  -  ygd  -  ^  -  y4(l  +  ^ 

=  >'1(1  +  I?)  -  y2<^  +  ?j)  -  VjCl  -  T))  +  y/l  -  rj) 

\.lo) 

The  integrals  are  evaluated  using  a  two -point  Gaussian  quadrature  in 
each  of  the  coordinates.  By  letting  k  =  1,2, 3, 4,  Eq.  (16)  gives  four  of 
the  equations  in  the  elemental  system.  The  other  eight  are  obtained 
from  Eqs.  (4a)  and  (4c). 

Artifical  viscosity  terms  were  included  at  times  in  order  to  main¬ 
tain  stability  in  the  transient  solution.  Since  the  magnitude  of  the  terms 
was  kept  small,  the  effect  on  the  final  steady-state  solution  is  assumed 
to  be  minor.  Thus,  it  was  assumed  that  an  accurate  evaluation  of  the 
viscous  terms  was  not  required  as  long  as  stability  was  maintained. 

The  following  describes  the  evaluation  of  the  viscous  terms  as  used 
in  this  report.  Consider  the  viscous  term  in  the  axial  momentum  equa¬ 
tion  without  the  constant  multiplier.  Using  Green's  theorem  results  in 

//  =  -  /f  n’'  .  Vudxdy  +  £  .  ^uds 


The  line  integral  was  ignored  giving 

-//j^  ■  VudA  =  _  i  rj  (nkfi'  +  J2’a'^)dxdy 

m  1=1 

Two-point  Gaussian  quadrature  was  also  used  to  approximate  this 
term.  The  quadrature  is  not  exact  in  this  case  because  of  the  Jacobian 
in  the  denominator, 

2.7  INTEGRATION  IN  TIME 

Equation  (11)  is  algebraically  linear  in  the  time  derivatives.  These 
derivatives  were  obtained  by  use  of  Cholesky  factorization. 
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Since  the  matrix  coefficients  are  not  functions  of  the  dependent  vari¬ 
ables,  the  factorization  need  be  performed  only  once.  The  solution  of 
the  system  of  equations  results  in  evaluation  of  the  time  derivatives 
which  can  be  used  in  a  numerical  integration  technique.  The  method 
used  by  the  authors  was  fourth-order-accurate  Runge-Kutta. 

Initial  conditions  were  obtained  from  a  potential  flow  program. 

The  program  solved  the  Cauchy-Riemann  equations  using  a  finite  ele¬ 
ment  method.  The  elements  used  were  the  same  in  both  programs.  In 
the  case  of  the  Cauchy-Riemann  equations,  the  resiQting  FEM  system 
is  a  linear  algebraic  system  which  can  be  solved  directly  without  iter¬ 
ation,  The  solution  to  the  Cauchy-Riemann  equations  gave  initial  con¬ 
ditions  for  u  and  v.  The  initial  conditions  for  T  were  obtained  by  using 
Bernoulli's  equation  in  the  following  form: 

T  =  I  -  _  1) 

\  2  /  (20) 

at  each  node  in  the  flow  field. 

Direct  numerical  integration  of  the  system  gave  unstable  results. 
From  observation  of  the  transient  solution,  it  was  noticed  that  the 
dependent  variable  T  was  the  first  to  go  unstable.  The  integration  pro¬ 
cedure  was  modified  to  reset  the  variable  T  after  a  specific  number  of 
integration  steps  had  been  taken.  The  u  and  v  variables  from  the  last 
integration  step  were  used  in  Eq.  (20)  to  obtain  the  T  variable  for  a  re¬ 
start.  This  modification  will  not  affect  the  converged  results  since 
Eq,  (20)  is  valid  at  steady  state.  For  the  subsonic  cases,  this  procedure 
resulted  in  a  stable  convergence  to  the  steady-state  solution.  The  best 
results  were  obtained  using  about  five  integration  steps  before  resetting 
T. 


In  the  cases  where  a  shock  was  present,  the  above  modifications 
were  not  sufficient.  For  these  cases  an  artificial  viscosity  term  was 
added  to  the  momentum  equation.  The  viscosity  term  was  kept  as  small 
as  possible  without  causing  an  instability.  Best  results  were  obtained 
by  systematically  reducing  the  amount  of  viscosity  as  the  solution 
started  to  converge.  In  these  cases  the  steady-state  conditions  were 
recognized  as  a  steady  cyclic  pattern  occurring  with  the  resetting  of 
the  initial  conditions.  This  behavior  is  caused  by  the  inconsistency 
between  the  viscous  term  in  the  momentum  equation  and  the  invlscid 
assumption  implied  in  Eq.  (20), 
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Although  only  the  banded  portion  of  the  linear  system  of  equations 
is  required,  the  storage  of  these  coefficients  can  be  an  obstacle  for 
large  grid  networks*  The  linear  system  was  modified  by  using  only 
the  center  threee  diagonals;  the  other  coefficients  are  set  equal  to  zero* 
The  resulting  tridiagonal  system  not  only  required  less  core  but 
resulted  in  faster  convergence  to  the  steady-state  solution.  This 
modification,  while  altering  the  transient  solution,  results  in  the  same 
steady-state  solution. 


3.0  RESULTS  AND  DISCUSSION 


Several  selected  problems  have  been  solved  in  order  to  verify  the 
numerical  scheme  discussed  in  previous  sections.  Subsonic  and  tran¬ 
sonic  flows  over  two-dimensional  airfoils  were  investigated.  The  air¬ 
foils  chosen  were  (1)  a  6-percent  biconvex  circular  arc,  (2)  an 
18-percent  biconvex  circular  arc,  and  (3)  the  NACA  nonlifting  four¬ 
digit  series  (NACA0012  and  NACA0024)  of  Abbott  et  al.  (Ref.  6). 

The  method  was  implemented  on  an  IBM  370/165  digital  computer. 
Storage  requirements  for  the  grid  networks  used  in  this  report  were 
155K  bytes  of  core.  A  heuristic  approach  was  taken  in  determining 
when  convergence  was  obtained.  Integration  was  continued  until  the  Cp 
along  the  center  streamline  held  steady  within  plot  accuracy  (approxi¬ 
mately  +  0.  001  to  +0.  01  depending  on  the  case).  For  subsonic  cases 
about  three  minutes  were  required.  When  a  shock  occurred,  twenty  to 
fifty  minutes  were  required,  depending  upon  the  strength  of  the  shock. 
Several  cases  were  run  for  additional  time  to  confirm  the  authors' 
judgment  in  determining  when  steady  state  occurs.  The  final  value 
of  the  artificial  viscosity,  a,  used  was  the  smallest  value  that  did  not 
give  wiggles  in  the  region  ahead  of  the  shock. 

The  predicted  pressure  distributions  for  various  free-stream 
Mach  numbers  have  been  correlated  with  experimental  data  and  with 
available  numerical  solutions.  Figure  3  shows  the  predicted  pressure 
distribution  along  the  axis  and  over  a  6 -percent  circular-arc  airfoil 
with  a  free-stream  Mach  number,  M^,  of  0.  7,  A  summary  of  pre¬ 
dicted  pressure  distributions,  compared  to  the  experimental  data  of 
Collins  and  Krupp  (Ref.  7)  and  Knechtel  (Ref,  8)  for  three  subcritical 
free-stream  Mach  numbers  (M„  =  0.710,  0.809,  and  0.  826)  and  a 
slightly  supercritical  case  (M„  =  0.  857),  is  shown  in  Fig.  4.  Further 
increase  in  free-stream  Mach  number  to  M^^  =  0.903  produced  a  super¬ 
critical  case  and  a  shock  as  shown  in  Fig.  5.  For  the  cases  where 
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M„  =  0.  857  and  0.  903,  results  obtained  by  Murman  (Refs.  9  and  10)  are 
also  shown.  For  the  lower  Mach  numbers  the  differences  between  the 
predictions  are  too  small  to  show  in  the  figure.  In  these  supercritical 
cases,  an  artificial  viscosity  (or  =  0.  005)  is  required  to  achieve  stability, 
which  has  an  effect  of  smearing  the  discontinuity  (Ref,  11).  The  pre¬ 
dicted  shock  location  being  upstream  of  Collin’s  data  is  attributed  to  the 
fixed  wall. 

Computed  Mach  number  contours  for  the  6 -percent  circular-arc 
airfoil  with  M^,  =  0. 903  are  shown  in  Fig,  6,  The  rough  contours  are 
attributed  to  poor  plotter  resolution  and  linear  interpolation  over  the 
sparse  grid  shown  in  Fig.  1, 

The  case  of  transonic  flow  over  an  18 -percent  circular-arc  airfoil 
for  =  0.741  is  shown  in  Fig.  7.  The  predicted  pressure  distribution 
over  the  airfoil  is  compared  to  the  experimental  data  of  McDevitt 
(Ref.  12),  A  very  weak  shock  is  predicted,  which  requires  only  a  very 
small  artificial  viscosity  (o  =  0.  001)  for  stability. 

Application  of  the  scheme  to  two  nonlifting  airfoils  of  the  NACA 
four-digit  series  {NACA0012  and  NACA0024,  Ref.  6)  is  shown  in 
Figs.  8,  9,  and  10,  Predictions  for  the  NACA0012  airfoil  were  calcu¬ 
lated  for  the  cases  of  M^  =  0.  72  and  =  0.  80.  Figures  8  and  9  show 
the  correlations  of  the  predicted  pressure  distributions  with  the  experi¬ 
mental  data  of  Vidal  et  al.  (Ref.  13).  Figure  2  shows  the  grid  used  for 
the  NACA0012  airfoil. 

A  low  Mach  number  case  for  the  NACA0024  airfoil  is  shown  in 
Fig.  10  for  the  comparison  of  predicted  velocity  distributions  with  the 
theoretical  calculations  of  Abbott  (Ref.  6). 

The  grid  network  used  for  these  problems  is  not  ideal  for  bl^lnt- 
leading-edge  configurations.  The  wiggles  present  in  the  NACA  wing 
solutions  are  stable  in  time  and  could  be  improved  with  a  more  refined 
grid  network. 


4.0  CONCLUDING  REMARKS 


The  finite  element  method  (FEM),  based  on  the  Galerkin  weighted 
residual  approach,  has  been  used  to  compute  flows  over  airfoils.  This 
method  was  found  to  be  sufficiently  accurate  and  reasonably  fast  for 
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subsonic  flow  over  circular-arc  airfoils.  However,  difficulty  was  en¬ 
countered  in  obtaining  stability  when  a  shock  occurred  in  supercritical 
cases.  In  the  present  approach,  the  problem  was  solved  through  the 
use  of  an  artifical  viscosity.  An  increase  in  computing  time  and  the 
shock-smearing  effect  were  observed  (Ref,  11).  Also,  the  NACA  air¬ 
foils  with  blunt  leading  edges  gave  difficulty  and  required  larger  con¬ 
vergence  times. 

This  is  a  preliminary  effort  on  the  application  of  the  FEM  to  the 
Euler  equations  in  the  transonic  flow  regime.  Additional  effort  will  be 
required  to  establish  this  scheme  as  an  engineering  tool  for  solutions 
of  fluid- dynamic  problems  of  the  type  considered.  The  convergence 
time  must  be  reduced  for  supercritical  cases  which  require  an  arti¬ 
ficial  viscosity.  Also,  in  future  work,  the  artificial  viscous  term 
should  be  replaced  with  the  correct  viscous  term  and  the  energy  equa¬ 
tion  should  be  added  to  the  system. 
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Figure  4.  Comparison  of  chordwise  pressure  distributions  for  a  6-percent 
biconvex  airfoil  as  a  function  of  Mach  number. 
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Figure  5.  Comparison  of  chordwise  pressure  distribution  for  supercritical  flow 
over  a  6-percent  biconvex  airfoil  at  =  0.903. 
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Axial  Distance  from  Leading  Edge,  x/C 

Figure  7.  Comparison  of  chordwise  pressure  distribution  over  an  18-percent 
biconvex  airfoil  at  »  0.741, 
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APPENDIX  A 

FINITE  ELEMENT  ALGORITHM  FOR  THE  EULERIAN  EQUATIONS 


Following  the  scheme  outlined  in  Section  2.0,  the  finite  element 
algorithm  is  developed  for  the  Eulerian  equations.  The  Galerkin 
method  of  weighted  residuals  is  used  to  obtain  a  system  of  ordinary 
differential  equations  in  time. 

Equation  (6)  is  substituted  into  Eq.  (5)  and  then  weighted  integrals 
are  formed  over  the  domain 
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where  k  =  1,  2,  3,  4 


Equation  (A-1)  can  be  rewritten  as 


-  1  STjU://  (O')  dxdy 

i=l i=l  ^  R_  * 


-  (y  -  1)  1  2  T;  V.  rj  n’^  dxdy 

i  =  n=i 

-  1  S  T.  V,  //  OMfl’)n‘'dxdy 
i=ii-i  ^  K  ^ 

(A-2a) 


2(f/  n^fi^dxd^d, 

i=i\  Rm  7  ■ 


♦  .  •  , 

2  2u.uJ/  n'^dxdy 


4  4 


2  2  V 

i=i)=i 


] 


Uj  //  dxdy 


4  .  , 

2  Tj  //  dxdy 

>  =  1  Rtn 


(A-2b) 


4  4 

2  2 
i=lj=l 


Vj//  dxdy 

'm 


4  4 

-  2  2viv:;/  dxdy 

i=ii=i  '  y 


(A-2c) 


The  above  is  a  system  of  twelve  equations  in  twelve  unknown  de¬ 
rivatives;  the  coefficients  are  evaluated  for  each  element  and  then 
assembled  to  form  Eq.  (11). 
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NOMENCLATURE 


C 
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CO 


N 

n 

P 

R 

S 

T 

t 

u 

•  •  • 

Ui,  Vi,  Ti 

V 

V 

W 

X 

y 

a 

y 

e 

?,r) 

p 

0 


Chord  Length 
Pressure  coefficient 
Jacobian 

7Pj((r  -  l)P„(u2)) 

Reference  Length 

Free -stream  Mach  Number 

Number  of  node  points 

Unit  normal  vector 

Pressure 

Domain 

Element  domain 
Arc  length 
Temperature 
Time 

Axial  component  of  velocity 
Time  derivatives  at  the  i^^  node 
Volume 

Transverse  component  of  velocity 
+  v^ 

Axial  coordinate  (measured  from  leading  edge) 

Transverse  coordinate 

Artificial  viscosity  coefficient 

Ratio  of  specific  heats 

Body  (airfoil)  angle 

Parametric  coordinates 

Density 

Basic  functions 
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Interpolation  function 
Del  operator 

Element  domain  boundary 


V 

8Rin 


SUBSCRIPTS 

t 

X 

y 

n 


Partial  derivative  with  respect  to  time 
Partial  derivative  with  respect  to  x 
Partial  derivative  with  respect  to  y 
Partial  derivative  with  respect  to  rj 
Partial  derivative  with  respect  to  | 
Free-stream  condition 


SUPERSCRIPTS 

*  Dimensional  parameter 
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